ENMTools
This package implements various tests, visualizations, and metrics for use with environmental niche models (ENMs) and species distribution models (SDMs).
Installation
At present, ENMTools is downloadable from https://github.com/danlwarren/ENMTools. There are multiple ways to download it. The easiest is to use devtools and install from GitHub.
Installing from GitHub using devtools
Run the following code from your R console:
install.packages("devtools")
library(devtools)
install_github("danlwarren/ENMTools")
library(ENMTools)Install from zip file
A zipped version of the package is available at https://github.com/danlwarren/ENMTools/archive/master.zip. To install from the zip file, download a copy of it to your system. Once it’s finished downloading, type the following (where PATH is the path to the zip file):
install.packages("devtools")
library(devtools)
install_local("PATH")
library(ENMTools)Interacting with ENMTools
Creating enmtools.species objects
First we’re going to load in some environmental data.
env.files <- list.files(path = "test/testdata/", pattern = "pc", full.names = TRUE)
env <- stack(env.files)
names(env) <- c("layer.1", "layer.2", "layer.3", "layer.4")
env <- setMinMax(env)ENMTools is primarily designed to examine patterns of similarity and difference between ENMs for different species. In order to simplify interactions with the functions in ENMTools, you need to put your data for each of your species into an enmtools.species object. You can create and view an empty enmtools.species object just by typing:
ahli <- enmtools.species()
ahli##
##
## Range raster not defined.
##
## Presence points not defined.
##
## Background points not defined.
##
## Species name not defined.You can add data to this object manually:
names(ahli)## [1] "range" "presence.points" "background.points"
## [4] "models" "species.name"ahli$species.name <- "ahli"
ahli$presence.points <- read.csv("test/testdata/ahli.csv")[,3:4]
ahli$range <- background.raster.buffer(ahli$presence.points, 50000, mask = env)
ahli$background.points <- background.points.buffer(points = ahli$presence.points,
radius = 20000, n = 1000, mask = env[[1]])
ahli##
##
## Range raster:
## class : RasterLayer
## dimensions : 418, 1535, 641630 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -86.90809, -74.11642, 19.80837, 23.2917 (xmin, xmax, ymin, ymax)
## coord. ref. : NA
## data source : in memory
## names : layer.1
## values : 1, 1 (min, max)
##
##
##
## Presence points (first ten only):
##
## Longitude Latitude
## ---------- ---------
## -80.0106 21.8744
## -79.9086 21.8095
## -79.8065 21.7631
## -79.8251 21.8095
## -79.8807 21.8374
## -79.9550 21.8374
## -80.3446 22.0136
## -80.2983 21.9951
## -80.1776 21.9023
## -80.1591 21.9673
##
##
## Background points (first ten only):
##
## Longitude Latitude
## ---------- ---------
## -80.43726 22.17087
## -79.99559 22.02920
## -79.89559 21.76254
## -79.76226 21.92920
## -80.22059 22.07087
## -80.36226 22.01254
## -80.01226 21.91254
## -80.48726 22.11254
## -80.07892 21.86254
## -80.32059 21.96254
##
##
## Species name: ahliOr you can add bits of it when the object is created:
allogus <- enmtools.species(species.name = "allogus",
presence.points = read.csv("test/testdata/allogus.csv")[,3:4])
allogus$range <- background.raster.buffer(allogus$presence.points, 50000, mask = env)
allogus$background.points <- background.points.buffer(points = allogus$presence.points,
radius = 20000, n = 1000, mask = env[[1]])
allogus##
##
## Range raster:
## class : RasterLayer
## dimensions : 418, 1535, 641630 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -86.90809, -74.11642, 19.80837, 23.2917 (xmin, xmax, ymin, ymax)
## coord. ref. : NA
## data source : in memory
## names : layer.1
## values : 1, 1 (min, max)
##
##
##
## Presence points (first ten only):
##
## Longitude Latitude
## ---------- ---------
## -79.2527 22.2109
## -78.7774 22.2241
## -78.6189 22.2373
## -78.1039 21.1809
## -78.0247 21.1809
## -77.9983 20.9301
## -77.9719 21.7091
## -77.9719 21.5507
## -77.9323 21.6167
## -77.9323 20.7320
##
##
## Background points (first ten only):
##
## Longitude Latitude
## ---------- ---------
## -75.28726 20.02920
## -75.86226 20.18754
## -76.66226 20.38754
## -74.38726 20.10420
## -77.82892 20.87087
## -78.08726 21.70420
## -76.32059 20.16254
## -76.63726 20.42087
## -76.97892 20.01254
## -78.88726 22.25420
##
##
## Species name: allogusBuilding an ENM
ENMTools contains functions to simplify the ENM construction process. Using enmtools.species objects and the correct modeling commands, we can build models very quickly. These commands are primarily wrappers to dismo model construction and projection functions, and at present are only available for GLM, Maxent, Domain, and Bioclim models. One of the nice bits about this setup is that it allows enmtools to automatically generate suitability maps, do model evaluation, and plot the marginal suitability of habitat for each variable separately.
GLM
GLMs usually require the user to supply a formula, an enmtools.species object, and some environmental data. If your formula is a strictly additive function of all of the environmental layers in env, though, enmtools.glm will build a formula automatically.
ahli.glm <- enmtools.glm(species = ahli, env = env, f = pres ~ layer.1 + layer.2 + layer.3 + layer.4, test.prop = 0.2)## Adding environmental data to species ahli
## Processing presence points...
## Processing background points...allogus.glm <- enmtools.glm(species = allogus, env = env, f = pres ~ layer.1 + layer.2 + layer.3 + layer.4, test.prop = 0.2)## Adding environmental data to species allogus
## Processing presence points...
## Processing background points...ahli.glm##
##
## Formula: presence ~ layer.1 + layer.2 + layer.3 + layer.4
## <environment: 0x7fc5b3dd05b8>
##
##
## Data table (top ten lines):
##
## Longitude Latitude layer.1 layer.2 layer.3 layer.4 presence
## --- ---------- --------- -------- -------- -------- -------- ---------
## 1 -80.0106 21.8744 2765 1235 1174 252 1
## 4 -79.8251 21.8095 2207 1877 967 259 1
## 5 -79.8807 21.8374 2244 1828 945 249 1
## 6 -79.9550 21.8374 2250 1766 919 235 1
## 7 -80.3446 22.0136 2201 1822 978 277 1
## 8 -80.2983 21.9951 2214 1786 986 284 1
## 9 -80.1776 21.9023 2287 1722 992 266 1
## 10 -80.1591 21.9673 2984 965 1311 237 1
## 11 -80.1498 21.9858 3042 841 1371 221 1
## 12 -80.1220 21.9301 2898 1033 1231 242 1
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)])
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.4711 -0.1761 -0.1204 -0.0845 3.2446
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 52.8490333 28.5760127 1.849 0.0644 .
## layer.1 -0.0144125 0.0070933 -2.032 0.0422 *
## layer.2 -0.0148595 0.0076379 -1.945 0.0517 .
## layer.3 -0.0006334 0.0080180 -0.079 0.9370
## layer.4 0.0085116 0.0248357 0.343 0.7318
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 130.29 on 1011 degrees of freedom
## Residual deviance: 119.42 on 1007 degrees of freedom
## AIC: 129.42
##
## Number of Fisher Scoring iterations: 8
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 12
## n absences : 1000
## AUC : 0.7642083
## cor : 0.09845587
## max TPR+TNR at : -4.40285
##
##
## Proportion of data wittheld for model fitting: 0.2
##
## Model fit (test data): class : ModelEvaluation
## n presences : 4
## n absences : 1000
## AUC : 0.743875
## cor : 0.05039809
## max TPR+TNR at : -4.517831
##
##
## Suitability:
## class : RasterLayer
## dimensions : 418, 1535, 641630 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -86.90809, -74.11642, 19.80837, 23.2917 (xmin, xmax, ymin, ymax)
## coord. ref. : NA
## data source : in memory
## names : layer
## values : 2.691814e-06, 0.9999737 (min, max)
Notice this produces the same formula as:
ahli.glm <- enmtools.glm(species = ahli, env = env, test.prop = 0.2)## Adding environmental data to species ahli
## Processing presence points...
## Processing background points...ahli.glm##
##
## Formula: presence ~ layer.1 + layer.2 + layer.3 + layer.4
## <environment: 0x7fc5b2c608e0>
##
##
## Data table (top ten lines):
##
## Longitude Latitude layer.1 layer.2 layer.3 layer.4 presence
## --- ---------- --------- -------- -------- -------- -------- ---------
## 1 -80.0106 21.8744 2765 1235 1174 252 1
## 2 -79.9086 21.8095 2289 1732 957 231 1
## 3 -79.8065 21.7631 2158 1870 983 253 1
## 5 -79.8807 21.8374 2244 1828 945 249 1
## 7 -80.3446 22.0136 2201 1822 978 277 1
## 8 -80.2983 21.9951 2214 1786 986 284 1
## 9 -80.1776 21.9023 2287 1722 992 266 1
## 10 -80.1591 21.9673 2984 965 1311 237 1
## 12 -80.1220 21.9301 2898 1033 1231 242 1
## 14 -80.2148 21.9394 2329 1692 1018 269 1
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)])
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.3660 -0.1837 -0.1240 -0.0843 3.1966
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 59.877906 27.899368 2.146 0.0319 *
## layer.1 -0.016190 0.006936 -2.334 0.0196 *
## layer.2 -0.016814 0.007426 -2.264 0.0236 *
## layer.3 -0.002510 0.008194 -0.306 0.7593
## layer.4 0.017902 0.023473 0.763 0.4457
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 130.29 on 1011 degrees of freedom
## Residual deviance: 120.92 on 1007 degrees of freedom
## AIC: 130.92
##
## Number of Fisher Scoring iterations: 8
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 12
## n absences : 1000
## AUC : 0.754
## cor : 0.08787457
## max TPR+TNR at : -4.56992
##
##
## Proportion of data wittheld for model fitting: 0.2
##
## Model fit (test data): class : ModelEvaluation
## n presences : 4
## n absences : 1000
## AUC : 0.74125
## cor : 0.05397096
## max TPR+TNR at : -3.360277
##
##
## Suitability:
## class : RasterLayer
## dimensions : 418, 1535, 641630 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -86.90809, -74.11642, 19.80837, 23.2917 (xmin, xmax, ymin, ymax)
## coord. ref. : NA
## data source : in memory
## names : layer
## values : 1.983417e-06, 0.9999955 (min, max)
If you want a more complicated formula, though (e.g., with interactions or polynomial effects), you’ll need to supply that manually.
ahli.glm <- enmtools.glm(species = ahli, env = env, f = pres ~ poly(layer.1, 2) + poly(layer.2, 2) + poly(layer.3, 2) + poly(layer.4, 2), test.prop = 0.2)## Adding environmental data to species ahli
## Processing presence points...
## Processing background points...ahli.glm##
##
## Formula: presence ~ poly(layer.1, 2) + poly(layer.2, 2) + poly(layer.3,
## 2) + poly(layer.4, 2)
## <environment: 0x7fc5ccf57238>
##
##
## Data table (top ten lines):
##
## Longitude Latitude layer.1 layer.2 layer.3 layer.4 presence
## --- ---------- --------- -------- -------- -------- -------- ---------
## 2 -79.9086 21.8095 2289 1732 957 231 1
## 3 -79.8065 21.7631 2158 1870 983 253 1
## 4 -79.8251 21.8095 2207 1877 967 259 1
## 6 -79.9550 21.8374 2250 1766 919 235 1
## 7 -80.3446 22.0136 2201 1822 978 277 1
## 8 -80.2983 21.9951 2214 1786 986 284 1
## 10 -80.1591 21.9673 2984 965 1311 237 1
## 11 -80.1498 21.9858 3042 841 1371 221 1
## 12 -80.1220 21.9301 2898 1033 1231 242 1
## 13 -80.1776 21.9673 2914 1020 1256 237 1
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)])
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.65758 -0.16717 -0.11277 -0.07515 3.09663
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.0894 0.4923 -10.337 <2e-16 ***
## poly(layer.1, 2)1 -111.7489 70.4889 -1.585 0.113
## poly(layer.1, 2)2 9.5637 28.0582 0.341 0.733
## poly(layer.2, 2)1 -114.3667 79.0120 -1.447 0.148
## poly(layer.2, 2)2 -5.6030 33.8869 -0.165 0.869
## poly(layer.3, 2)1 2.0866 44.4962 0.047 0.963
## poly(layer.3, 2)2 11.8612 17.7449 0.668 0.504
## poly(layer.4, 2)1 3.3925 15.1844 0.223 0.823
## poly(layer.4, 2)2 -7.8530 12.9228 -0.608 0.543
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 130.29 on 1011 degrees of freedom
## Residual deviance: 114.97 on 1003 degrees of freedom
## AIC: 132.97
##
## Number of Fisher Scoring iterations: 8
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 12
## n absences : 1000
## AUC : 0.8186667
## cor : 0.118871
## max TPR+TNR at : -4.786353
##
##
## Proportion of data wittheld for model fitting: 0.2
##
## Model fit (test data): class : ModelEvaluation
## n presences : 4
## n absences : 1000
## AUC : 0.68675
## cor : 0.03338417
## max TPR+TNR at : -4.718923
##
##
## Suitability:
## class : RasterLayer
## dimensions : 418, 1535, 641630 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -86.90809, -74.11642, 19.80837, 23.2917 (xmin, xmax, ymin, ymax)
## coord. ref. : NA
## data source : in memory
## names : layer
## values : 2.220446e-16, 0.9999816 (min, max)
To check out the marginal response functions, you only need to type
ahli.glm$response.plots## $layer.1
##
## $layer.2
##
## $layer.3
##
## $layer.4
You can also visualize your models and data in a 2D environment space using any pair of layers from your environment stack. These plots hold all non-plotted variables (layer.1 and layer.3 in this case) constant at their mean value across all presence points, then vary the plotted variables between the minimum and maximum values in env.
The suit.plot shows you suitability in environment space as a function of your two variables, with brighter colors representing variable combinations predicted to be more suitable. The points represent the occurrence points for your species in that environment space.
The colored raster of the background.plot shows you the density of background points in environment space, while the white points again represent your occurrence points in environment space.
visualize.enm(ahli.glm, env, layers = c("layer.2", "layer.4"))## $suit.plot
##
## $background.plot
GAM, Bioclim, Domain, and Maxent
The procedure for building Bioclim, Domain, and Maxent models is similar to the procedure for GLMs, with the exception that you do not need to pass a formula to the model function. Note that running Maxent models requires a bit of extra setup; see dismo documentation for details.
ahli.gam <- enmtools.gam(ahli, env, test.prop = 0.2)
ahli.dm <- enmtools.dm(ahli, env, test.prop = 0.2)
ahli.bc <- enmtools.bc(ahli, env, test.prop = 0.2)
ahli.mx <- enmtools.maxent(ahli, env, test.prop = 0.2)Metrics: breadth, correlation, and overlap
ENMTools provides a number of metrics for ENMs and for similarities between ENMs. These include measures of niche breadth, based on Levins(1968). An important caveat when interpreting these metrics is that they are driven to some (variable) extent by the availability of different combinations of environmental variables. As such they are more accurately interpreted as a measurment of the smoothness of the geographic distribution of suitability scores than as an estimate of the breadth of the fundamental niche; an organism with narrow fundamental niche breadth that nonetheless encompasses a set of environmental conditions that is quite common will have a high breadth when measured using ENMs, while having a low breadth in environment space.
raster.breadth(ahli.glm)## $B1
## [1] 0.8650396
##
## $B2
## [1] 0.1802174ENMTools also provides metrics for measuring similarity between ENMs. These include Schoener’s D (Schoener 1968), I (Warren et al. 2008), and the Spearman rank correlation coefficient between two rasters. While D and I are commonly used in the ENM literature, they may tend to overestimate similarity between ENMs when many grid cells are of similar values (e.g., when two species prefer different habitat but the region contains a great deal of habitat that is unsuitable for both).
raster.overlap(ahli.glm, allogus.glm)## $D
## [1] 0.3616807
##
## $I
## [1] 0.6267567
##
## $rank.cor
## [1] 0.4740915A new feature of the R version of ENMTools is that you can now use these same metrics in the n-dimensional space of all combinations of environmental variables, instead of restricting your measures of model similarity to those sets of conditions that appear in the training region. This is done by repeatedly drawing Latin hypercube samples from the space of all possible combinations of environmental variables given the min and max of each variable within the training region. ENMTools continues to draw samples until subsequent iterations differ by less than a specified tolerance value. Lower tolerance values result in more precise estimates of overlap, but can take much longer to calculate.
env.overlap(ahli.glm, allogus.glm, env, tolerance = .001)## $env.D
## [1] 0.2753451
##
## $env.I
## [1] 0.4720855
##
## $env.cor
## [1] -0.00658996Hypothesis testing
Niche identity or equivalency test
In this example, we will run a niche identity (also called equivalency) test, as in Warren et al. 2008. This test takes the presence points for a pair of species and randomly reassigns them to each species, then builds ENMs for these randomized occurrences. By doing this many times, we can estimate the probability distribution for ENM overlap between species under the null hypothesis that the two species’ occurrences in the environment are effectively a random draw from the same underlying distribution. Note that niche evolution is only one of many reasons why two species’ realized environmental distributions might cause departures from this null hypothesis. See Warren et al. 2014 for details.
To run an identity test, we need to decide what type of models we will build, how many replicates we will run, and (in the case of GLM) a model formula to use for empirical models and the Monte Carlo replicates. The resulting object contains the replicate models, p values, and plots of the results. Typically idenity tests are run with at least 99 replicates, but we are using a smaller number here for the sake of execution time.
NOTE: In order for models to be comparable, both empirical and pseudoreplicate models for the identity test are conducted with pseudoabsence points pooled for the two species being compared.
id.glm <- identity.test(species.1 = ahli, species.2 = allogus, env = env, type = "glm", nreps = 4)id.glm##
##
##
##
## Identity test ahli vs. allogus
##
## Identity test p-values:
## D I rank.cor env.D env.I env.cor
## 0.2 0.2 0.2 0.2 0.2 0.2
##
##
## Replicates:
##
##
## D I rank.cor env.D env.I env.cor
## ---------- ---------- ---------- ----------- ---------- ---------- -----------
## empirical 0.2548982 0.5075499 -0.4020943 0.0050492 0.0299701 -0.5875723
## rep 1 0.8546444 0.9844519 0.7852184 0.7513282 0.9482161 0.8399446
## rep 2 0.7665889 0.9582296 0.3574021 0.6183091 0.8775111 0.6378208
## rep 3 0.8195256 0.9752231 0.9737111 0.7759048 0.9546163 0.9225750
## rep 4 0.9490372 0.9976710 0.9807814 0.9115978 0.9939383 0.9876159
Background or similarity test
The background or similarity test compares the overlap seen between two species’ ENMs to the overlap expected by chance if one or both species was effectively choosing habitat at random from within their broad geographic range. The purpose of this test is to correct for the availability of habitat and ask whether the observed similarity between species or populations is significantly more (or less) than expected given the available set of environments in the regions in which they occur.
NOTE: In order for models to be comparable, both empirical and pseudoreplicate models for the background test are conducted with pseudoabsence points pooled for the two species being compared.
In Warren et al. 2008, we developed this test in the context of comparing one species’ actual occurrence to the random background occurrences of the other species. This is what we call an “asymmetric” test, and in our case we did the test in both directions with the idea that we might compare the results of A vs. B background to the results of B vs. A background. This may be informative in some cases, but many people have also found this asymmetry confusing (and indeed it is often difficult to interpret). For that reason, the background test here can be conducted against a null hypothesis that is generated from “asymmetric” (species.1 vs species.2 background) or “symmetric” (species.1 background vs. species.2 background) comparisons.
Here, for instance, is a Bioclim background test using the classical asymmetric approach:
bg.bc.asym <- background.test(species.1 = ahli, species.2 = allogus, env = env, type = "bc", nreps = 4, test.type = "asymmetric" )bg.bc.asym##
##
##
##
## Asymmetric background test ahli vs. allogus background
##
## background test p-values:
## D I rank.cor env.D env.I env.cor
## 0.2 0.2 0.2 0.2 0.2 0.2
##
##
## Replicates:
##
##
## D I rank.cor env.D env.I env.cor
## ---------- ---------- ---------- ---------- ---------- ---------- ----------
## empirical 0.1328502 0.3177390 0.0706201 0.0205455 0.1116701 0.0861670
## rep 1 0.1486335 0.3305719 0.1450342 0.0424154 0.1589587 0.0986768
## rep 2 0.1651580 0.3500385 0.1639466 0.0501244 0.1753446 0.1255899
## rep 3 0.1665392 0.3649020 0.1905156 0.0751387 0.2240976 0.1826939
## rep 4 0.1624940 0.3544069 0.1801364 0.0498466 0.1885051 0.1350149
And here is a Domain background test using the symmetric approach:
bg.dm.sym <- background.test(species.1 = ahli, species.2 = allogus, env = env, type = "dm", nreps = 4, test.type = "symmetric" )bg.dm.sym##
##
##
##
## Symmetric background test ahli background vs. allogus background
##
## background test p-values:
## D I rank.cor env.D env.I env.cor
## 0.2 0.2 0.2 0.2 0.2 0.2
##
##
## Replicates:
##
##
## D I rank.cor env.D env.I env.cor
## ---------- ---------- ---------- ---------- ---------- ---------- ----------
## empirical 0.4929334 0.7052122 0.2916150 0.1067514 0.3106242 0.2245317
## rep 1 0.9633144 0.9983875 0.8452877 0.7789384 0.9084278 0.8025493
## rep 2 0.9683391 0.9979490 0.9703954 0.5167449 0.7288396 0.7550513
## rep 3 0.9367693 0.9887789 0.8515679 0.2534965 0.4923326 0.4736000
## rep 4 0.9212382 0.9927110 0.6998317 0.6713370 0.8399510 0.7405201
Ecospat tests
Using enmtools.species objects also provides a simplified interface to the niche equivalency and similarity tests (or identity and background tests, respectively) that were developed by Broennimann et al. (2012). These tests do not rely on ENMs, instead using kernel density smoothing to estimate density of the species in environment space. Ecospat also uses the density of the available environment to correct for availability when measuring overlaps, so that overlaps are not strictly driven by availability of combinations of environmental variables.
These tests only work with two environmental axes, so they are often done with the top two PC axes of a set of environments. In our case we’ll just pick a couple of environmental layers, though (layer.1 and layer.2). Here’s an equivalency/identity test:
esp.id <- enmtools.ecospat.id(ahli, allogus, env[[c("layer.1", "layer.2")]])## ===========================================================================esp.id##
##
##
##
## Ecospat identity test ahli vs. allogus
##
## Species Longitude Latitude layer.1 layer.2
## -------- ---------- --------- -------- --------
## ahli -80.0106 21.8744 2765 1235
## ahli -79.9086 21.8095 2289 1732
## ahli -79.8065 21.7631 2158 1870
## ahli -79.8251 21.8095 2207 1877
## ahli -79.8807 21.8374 2244 1828
## ahli -79.9550 21.8374 2250 1766
##
##
## Species Longitude Latitude layer.1 layer.2
## -------- ---------- --------- -------- --------
## ahli.bg -80.43726 22.17087 2324 1828
## ahli.bg -79.99559 22.02920 2660 1382
## ahli.bg -79.89559 21.76254 2184 1849
## ahli.bg -79.76226 21.92920 2440 1733
## ahli.bg -80.22059 22.07087 2458 1589
## ahli.bg -80.36226 22.01254 2240 1817
##
##
## Species Longitude Latitude layer.1 layer.2
## ----------- ---------- --------- -------- --------
## allogus.bg -79.2527 22.2109 2656 1683
## allogus.bg -78.7774 22.2241 2373 1967
## allogus.bg -78.6189 22.2373 2317 1980
## allogus.bg -78.1039 21.1809 2461 1669
## allogus.bg -78.0247 21.1809 2437 1684
## allogus.bg -77.9983 20.9301 2315 1746
##
##
## Species Longitude Latitude layer.1 layer.2
## -------- ---------- --------- -------- --------
## allogus -75.28726 20.02920 1692 1774
## allogus -75.86226 20.18754 2181 1479
## allogus -76.66226 20.38754 2107 1720
## allogus -74.38726 20.10420 2046 1422
## allogus -77.82892 20.87087 2338 1703
## allogus -78.08726 21.70420 2380 1821
##
##
## Species Longitude Latitude layer.1 layer.2
## ----------- ---------- --------- -------- --------
## background -80.91226 23.27087 2495 1993
## background -80.90392 23.27087 2504 2000
## background -80.93726 23.26254 2500 1993
## background -80.92892 23.26254 2500 1993
## background -80.92059 23.26254 2509 2002
## background -80.91226 23.26254 2504 2000
##
##
## ecospat.id test p-values:
## D I
## 0.02 0.02
## TableGrob (3 x 2) "arrange": 6 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
## 5 5 (3-3,1-1) arrange gtable[layout]
## 6 6 (3-3,2-2) arrange gtable[layout]And here’s a symmetric background test. The difference between symmetric and asymmetric for these tests is the same as for the background tests presented above.
esp.bg.sym <- enmtools.ecospat.bg(ahli, allogus, env[[c("layer.1", "layer.3")]], test.type = "symmetric")## ===========================================================================esp.bg.sym##
##
##
##
## Ecospat background test symmetric ahli vs. allogus
##
## Species Longitude Latitude layer.1 layer.3
## -------- ---------- --------- -------- --------
## ahli -80.0106 21.8744 2765 1174
## ahli -79.9086 21.8095 2289 957
## ahli -79.8065 21.7631 2158 983
## ahli -79.8251 21.8095 2207 967
## ahli -79.8807 21.8374 2244 945
## ahli -79.9550 21.8374 2250 919
##
##
## Species Longitude Latitude layer.1 layer.3
## -------- ---------- --------- -------- --------
## ahli.bg -80.43726 22.17087 2324 938
## ahli.bg -79.99559 22.02920 2660 1077
## ahli.bg -79.89559 21.76254 2184 948
## ahli.bg -79.76226 21.92920 2440 1000
## ahli.bg -80.22059 22.07087 2458 1028
## ahli.bg -80.36226 22.01254 2240 1006
##
##
## Species Longitude Latitude layer.1 layer.3
## ----------- ---------- --------- -------- --------
## allogus.bg -79.2527 22.2109 2656 1097
## allogus.bg -78.7774 22.2241 2373 1067
## allogus.bg -78.6189 22.2373 2317 1065
## allogus.bg -78.1039 21.1809 2461 846
## allogus.bg -78.0247 21.1809 2437 877
## allogus.bg -77.9983 20.9301 2315 907
##
##
## Species Longitude Latitude layer.1 layer.3
## -------- ---------- --------- -------- --------
## allogus -75.28726 20.02920 1692 951
## allogus -75.86226 20.18754 2181 915
## allogus -76.66226 20.38754 2107 907
## allogus -74.38726 20.10420 2046 797
## allogus -77.82892 20.87087 2338 866
## allogus -78.08726 21.70420 2380 1002
##
##
## Species Longitude Latitude layer.1 layer.3
## ----------- ---------- --------- -------- --------
## background -80.91226 23.27087 2495 1075
## background -80.90392 23.27087 2504 1082
## background -80.93726 23.26254 2500 1077
## background -80.92892 23.26254 2500 1077
## background -80.92059 23.26254 2509 1082
## background -80.91226 23.26254 2504 1080
##
##
## ecospat.bg test p-values:
## D I
## 0.56 0.44
## TableGrob (3 x 2) "arrange": 6 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
## 5 5 (3-3,1-1) arrange gtable[layout]
## 6 6 (3-3,2-2) arrange gtable[layout]Rangebreak tests
ENMTools also allows you to perform linear, blob, and ribbon rangebreak tests as developed in Glor and Warren 2011. The linear and blob tests are two versions of a test that permit one to ask whether the geographic regions occupied by two species are more environmentally different than expected by chance. The ribbon test, meanwhile, is designed to test whether the ranges of two species are divided by a region that is relatively unsuitable to one or both forms.
For the linear and blob tests, you call them very much like you would the identity and background tests. Here’s a linear one using GLM models:
rbl.glm <- rangebreak.linear(ahli, allogus, env, type = "glm", nreps = 4)##
## Building empirical models...
## Adding environmental data to species ahli
## Processing presence points...
## Processing background points...## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Adding environmental data to species allogus
## Processing presence points...
## Processing background points...
##
## Building replicate models...
##
## Replicate 1 ...
## Adding environmental data to species ahli
## Processing presence points...
## Processing background points...## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Adding environmental data to species allogus
## Processing presence points...
## Processing background points...
##
## Replicate 2 ...
## Adding environmental data to species ahli
## Processing presence points...
## Processing background points...## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Adding environmental data to species allogus
## Processing presence points...
## Processing background points...
##
## Replicate 3 ...
## Adding environmental data to species ahli
## Processing presence points...
## Processing background points...
## Adding environmental data to species allogus
## Processing presence points...
## Processing background points...
##
## Replicate 4 ...
## Adding environmental data to species ahli
## Processing presence points...
## Processing background points...## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Adding environmental data to species allogus
## Processing presence points...
## Processing background points...rbl.glm##
##
##
##
## Linear rangebreak test ahli vs. allogus
##
## rangebreak test p-values:
## D I rank.cor env.D env.I env.cor
## 1.0 1.0 0.6 0.4 0.4 0.2
##
##
## Replicates:
##
##
## D I rank.cor env.D env.I env.cor
## ---------- ---------- ---------- ----------- ---------- ---------- -----------
## empirical 0.2548982 0.5075499 -0.4020943 0.0048436 0.0298245 -0.5920514
## rep 1 0.0726801 0.2408095 0.5951234 0.5345011 0.7621649 0.5297946
## rep 2 0.2548982 0.5075499 -0.4020943 0.0046822 0.0293862 -0.5909788
## rep 3 0.1811089 0.4455980 -0.1197877 0.3190094 0.5509033 -0.0320241
## rep 4 0.2548982 0.5075499 -0.4020943 0.0050295 0.0301655 -0.5847006
And here’s a blob test using Bioclim:
rbb.bc <- rangebreak.blob(ahli, allogus, env, type = "bc", nreps = 4)##
## Building empirical models...
##
## Building replicate models...
##
## Replicate 1 ...
##
## Replicate 2 ...
##
## Replicate 3 ...
##
## Replicate 4 ...rbb.bc##
##
##
##
## blob rangebreak test ahli vs. allogus
##
## rangebreak test p-values:
## D I rank.cor env.D env.I env.cor
## 0.8 0.8 0.8 0.4 0.4 0.4
##
##
## Replicates:
##
##
## D I rank.cor env.D env.I env.cor
## ---------- ---------- ---------- ----------- ---------- ---------- ----------
## empirical 0.1328502 0.3177390 0.0706201 0.0223614 0.1171762 0.0934946
## rep 1 0.0437390 0.1507199 -0.0311702 0.1773347 0.2482751 0.2017040
## rep 2 0.1328502 0.3177390 0.0706201 0.0198827 0.1070997 0.0823804
## rep 3 0.0327940 0.1544608 -0.0172882 0.2056749 0.2907126 0.2460994
## rep 4 0.1608184 0.3470140 0.1160445 0.3570956 0.5817798 0.5617702
If you want to access the individual replicates (for instance to see how your ranges are being split up), you can find them in the list named “replicate.models” inside your rangebreak test object.
rbl.glm$replicate.models$ahli.rep.1##
##
## Formula: presence ~ layer.1 + layer.2 + layer.3 + layer.4
## <environment: 0x7fc5b38076d0>
##
##
## Data table (top ten lines):
##
## Longitude Latitude layer.1 layer.2 layer.3 layer.4 presence
## --- ---------- --------- -------- -------- -------- -------- ---------
## 79 -74.2483 20.2963 1935 1503 923 521 1
## 64 -74.8161 20.5867 2558 1292 660 706 1
## 68 -74.6312 20.4679 2934 960 332 815 1
## 59 -74.9877 20.6528 2328 1492 779 664 1
## 61 -74.9481 20.6264 2345 1410 738 656 1
## 81 -74.1691 20.2170 1608 1722 972 445 1
## 50 -75.6611 20.9961 1828 1777 1041 497 1
## 71 -74.5916 20.4283 2944 930 364 805 1
## 78 -74.2483 20.2434 1933 1431 954 485 1
## 63 -74.8293 20.5471 2418 1234 903 604 1
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)])
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.97625 -0.00030 -0.00002 -0.00001 2.74322
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 21.915940 10.371697 2.113 0.03460 *
## layer.1 -0.017237 0.006085 -2.833 0.00462 **
## layer.2 -0.012974 0.005079 -2.554 0.01064 *
## layer.3 0.005104 0.004761 1.072 0.28377
## layer.4 0.047478 0.018312 2.593 0.00952 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 186.634 on 2015 degrees of freedom
## Residual deviance: 93.555 on 2011 degrees of freedom
## AIC: 103.55
##
## Number of Fisher Scoring iterations: 13
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 16
## n absences : 2000
## AUC : 0.980375
## cor : 0.2265177
## max TPR+TNR at : -3.739233
##
##
## Proportion of data wittheld for model fitting: 0
##
## Model fit (test data): [1] NA
##
##
## Suitability:
## class : RasterLayer
## dimensions : 418, 1535, 641630 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -86.90809, -74.11642, 19.80837, 23.2917 (xmin, xmax, ymin, ymax)
## coord. ref. : NA
## data source : in memory
## names : layer
## values : 2.220446e-16, 0.6496604 (min, max)
rbl.glm$replicate.models$allogus.rep.1##
##
## Formula: presence ~ layer.1 + layer.2 + layer.3 + layer.4
## <environment: 0x7fc5c82241c0>
##
##
## Data table (top ten lines):
##
## Longitude Latitude layer.1 layer.2 layer.3 layer.4 presence
## --- ---------- --------- -------- -------- -------- -------- ---------
## 62 -74.9349 20.5339 2312 1107 1243 427 1
## 69 -74.6308 20.3725 2497 1110 792 616 1
## 54 -75.2782 20.6396 2089 1642 871 617 1
## 76 -74.2747 20.1114 1972 1243 1112 422 1
## 72 -74.4595 20.1906 2190 1196 994 502 1
## 74 -74.3803 20.1378 2060 1291 918 503 1
## 65 -74.8029 20.3359 2268 1102 1199 419 1
## 19 -78.6189 22.2373 2317 1980 1065 374 1
## 70 -74.6180 20.1378 2173 1157 1082 452 1
## 51 -75.6611 20.6792 2207 1706 875 487 1
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)])
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.3716 -0.2328 -0.1894 -0.1726 3.0655
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.1017704 2.4440531 0.860 0.3898
## layer.1 -0.0011967 0.0005463 -2.190 0.0285 *
## layer.2 -0.0015443 0.0005335 -2.894 0.0038 **
## layer.3 -0.0005267 0.0012190 -0.432 0.6657
## layer.4 -0.0009379 0.0013137 -0.714 0.4753
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 629.56 on 3064 degrees of freedom
## Residual deviance: 618.18 on 3060 degrees of freedom
## AIC: 628.18
##
## Number of Fisher Scoring iterations: 7
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 65
## n absences : 3000
## AUC : 0.6338564
## cor : 0.06187162
## max TPR+TNR at : -3.875799
##
##
## Proportion of data wittheld for model fitting: 0
##
## Model fit (test data): [1] NA
##
##
## Suitability:
## class : RasterLayer
## dimensions : 418, 1535, 641630 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -86.90809, -74.11642, 19.80837, 23.2917 (xmin, xmax, ymin, ymax)
## coord. ref. : NA
## data source : in memory
## names : layer
## values : 0.006466657, 0.07945109 (min, max)
For the ribbon rangebreak test, you will need one extra thing; a third enmtools.species object representing the occurrence points (for one or both species) that fall within the ribbon of putatively unsuitable habitat. In the case of these two anoles we don’t have such a ribbon, so we’ll just simulate one based on some random points.
ribbon <- enmtools.species(species.name = "ribbon")
ribbon$presence.points <- data.frame(Longitude = runif(n = 10, min = -79, max = -78.5),
Latitude = runif(n = 10, min = 21.7, max = 22.1))
plot(env[[1]])
points(ribbon$presence.points, pch = 16)
ribbon$range <- background.raster.buffer(ribbon$presence.points, 20000, mask = env)
ribbon##
##
## Range raster:
## class : RasterLayer
## dimensions : 418, 1535, 641630 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -86.90809, -74.11642, 19.80837, 23.2917 (xmin, xmax, ymin, ymax)
## coord. ref. : NA
## data source : in memory
## names : layer.1
## values : 1, 1 (min, max)
##
##
##
## Presence points (first ten only):
##
## Longitude Latitude
## ---------- ---------
## -78.64187 21.85732
## -78.52183 21.79001
## -78.61741 21.71052
## -78.68714 21.77541
## -78.73830 22.01088
## -78.60156 22.07664
## -78.56221 21.87031
## -78.88530 21.85005
## -78.89807 21.90540
## -78.90946 21.92327
##
##
## Background points not defined.
##
## Species name: ribbonNow we’ll run a ribbon rangebreak test using GLM models with quadratic effects. We also need to tell it the width of the ribbons to generate for the replicates. The units for the width argument are the same units that the presence points are in; e.g., if the points are in decimal degrees you should supply the width of the barrier in decimal degrees.
rbr.glm <- rangebreak.ribbon(ahli, allogus, ribbon, env, type = "glm", f = pres ~ poly(layer.1, 2) + poly(layer.2, 2) + poly(layer.3, 2) + poly(layer.4, 2), width = 0.5, nreps = 4)##
##
## No background points provided, drawing background from range raster.## Warning in couldBeLonLat(mask): CRS is NA. Assuming it is longitude/
## latitude##
## Building empirical models...
## Adding environmental data to species ahli
## Processing presence points...
## Processing background points...## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Adding environmental data to species allogus
## Processing presence points...
## Processing background points...
## Adding environmental data to species ribbon
## Processing presence points...
## Processing background points...## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Adding environmental data to species outside
## Processing presence points...
## Processing background points...
##
## Building replicate models...
##
## Replicate 1 ...
## Adding environmental data to species ahli
## Processing presence points...
## Processing background points...## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Adding environmental data to species allogus
## Processing presence points...
## Processing background points...
## Adding environmental data to species ribbon
## Processing presence points...
## Processing background points...## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Adding environmental data to species outside
## Processing presence points...
## Processing background points...
##
## Replicate 2 ...
## Adding environmental data to species ahli
## Processing presence points...
## Processing background points...## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Adding environmental data to species allogus
## Processing presence points...
## Processing background points...
## Adding environmental data to species ribbon
## Processing presence points...
## Processing background points...
## Adding environmental data to species outside
## Processing presence points...
## Processing background points...
##
## Replicate 3 ...
## Adding environmental data to species ahli
## Processing presence points...
## Processing background points...## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Adding environmental data to species allogus
## Processing presence points...
## Processing background points...
## Adding environmental data to species ribbon
## Processing presence points...
## Processing background points...## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Adding environmental data to species outside
## Processing presence points...
## Processing background points...
##
## Replicate 4 ...
## Adding environmental data to species ahli
## Processing presence points...
## Processing background points...
## Adding environmental data to species allogus
## Processing presence points...
## Processing background points...
## Adding environmental data to species ribbon
## Processing presence points...
## Processing background points...## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Adding environmental data to species outside
## Processing presence points...
## Processing background points...rbr.glm##
##
##
##
## ribbon rangebreak test ahli vs. allogus
##
## rangebreak test p-values...
##
## Species 1 vs. Species 2:
## D I rank.cor env.D env.I env.cor
## 0.6 0.4 1.0 0.8 0.8 1.0
##
## Species 1 vs. Ribbon:
## D I rank.cor env.D env.I env.cor
## 0.4 0.4 0.2 0.4 0.4 0.4
##
## Species 2 vs. Ribbon:
## D I rank.cor env.D env.I env.cor
## 0.2 0.2 0.2 0.2 0.2 0.8
##
## Outside vs. Ribbon:
## D I rank.cor env.D env.I env.cor
## 0.4 0.4 0.2 0.2 0.2 0.8
##
##
## Replicates:
##
## Species 1 vs. Species 2:
## D I rank.cor env.D env.I
## empirical 0.2742489 0.5051508 0.15455193 0.020467893 0.07339637
## rep 1 0.2723466 0.5065487 0.07376484 0.011260859 0.05189974
## rep 2 0.1141439 0.3173775 -0.01915522 0.002053111 0.01687757
## rep 3 0.2785799 0.5187816 0.01615524 0.013228309 0.06217732
## rep 4 0.4299572 0.7161948 -0.18828103 0.149242318 0.30937170
## env.cor
## empirical -0.2556596
## rep 1 -0.4034983
## rep 2 -0.8029448
## rep 3 -0.4230139
## rep 4 -0.3242595
##
## Species 1 vs. Ribbon:
## D I rank.cor env.D env.I
## empirical 0.12405130 0.2641380 0.1828724 0.0007867978 0.006713102
## rep 1 0.33688792 0.5881668 0.3995109 0.0099571507 0.041027548
## rep 2 0.22134508 0.4768959 0.5665779 0.4399150549 0.666858437
## rep 3 0.32758922 0.5735703 0.3837598 0.0257900927 0.070656140
## rep 4 0.05345405 0.1904684 0.3416105 0.0001623207 0.005240783
## env.cor
## empirical 0.1477843
## rep 1 0.4589192
## rep 2 0.6446700
## rep 3 0.4641404
## rep 4 -0.2352111
##
## Species 2 vs. Ribbon:
## D I rank.cor env.D env.I
## empirical 0.05153580 0.1803862 -0.71437687 9.308834e-05 0.003681101
## rep 1 0.45363790 0.7438772 0.27883255 5.485034e-02 0.176110109
## rep 2 0.43150795 0.7351777 -0.06515655 3.705948e-02 0.118615914
## rep 3 0.46769553 0.7511526 0.34464961 4.797471e-02 0.163117177
## rep 4 0.06043074 0.1964880 0.14811053 3.090985e-02 0.147245259
## env.cor
## empirical -0.125132471
## rep 1 -0.280372937
## rep 2 -0.668879084
## rep 3 -0.256833051
## rep 4 -0.002399887
##
## Outside vs. Ribbon:
## D I rank.cor env.D env.I
## empirical 0.06687670 0.2131445 -0.73482153 0.0001481916 0.004870576
## rep 1 0.48150293 0.7682087 0.21838188 0.0453866971 0.149897383
## rep 2 0.56320319 0.8460272 0.66302443 0.1047377506 0.245976321
## rep 3 0.48355076 0.7643736 0.29185424 0.0384509537 0.133198826
## rep 4 0.06067813 0.1981604 0.02584452 0.0350447558 0.158986867
## env.cor
## empirical -0.140034
## rep 1 -0.472773
## rep 2 -0.420043
## rep 3 -0.451324
## rep 4 -0.069941
Note that the output table here has slope, intercept, and intercept offset.
rbr.glm$lines.df## slope intercept offset
## 1 -0.9582584 -53.211664 0.3462531
## 2 -0.2418844 2.152673 0.2572096
## 3 -1.1035993 -64.467252 0.3723180
## 4 0.9509825 94.811239 0.3449971The intercept denotes the intercept corresponding to the CENTER of each ribbon. To get the lines denoting the edges of the ribbons (for example if you want to plot the ribbons on a map), you add and substract the offset. In other words, the top edge of the ribbon is given by y = (slope * x) + intercept + offset, while the bottom edge is given by y = (slope * x) + intercept - offset.
Literature cited
Broennimann, O., Fitzpatrick, M. C., Pearman, P. B., Petitpierre, B., Pellissier, L., Yoccoz, N. G., Thuiller, W., Fortin, M.-J., Randin, C., Zimmermann, N. E., Graham, C. H. and Guisan, A. (2012), Measuring ecological niche overlap from occurrence and spatial environmental data. Global Ecology and Biogeography, 21: 481–497. doi:10.1111/j.1466-8238.2011.00698.x
Levins, R. 1968. Evolution In Changing Environments. Monographs in Population Biology, volume 2. Princeton University Press, Princeton, New Jersey, USA.
Schoener, T. W. 1968. Anolis lizards of Bimini: resource partitioning in a complex fauna. Ecology 49:704- 726.
Warren, D.L., R.E. Glor, and M. Turelli. 2008. Environmental niche identity versus conservatism: quantitative approaches to niche evolution. Evolution 62:2868-2883. doi: 10.1111/j.1558-5646.2008.00482.x
Warren, D.L., M. Cardillo, D.F. Rosauer, and D.I. Bolnick. 2014. Mistaking geography for biology: inferring processes from species distributions. Trends in Ecology and Evolution 29 (10), 572-580. doi: 10.1016/j.tree.2014.08.003